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Abstract.—Population genetic analyses are a powerful tool for obtaining information about cryptic genetic 
lineages, population structure, and the distribution of intra- and interpopulation genetic diversity across the 
landscape. This knowledge is crucial for establishing units for the conservation management of endangered 
species. Species with limited dispersal capacities, such as amphibians, are particularly affected by habitat 
fragmentation and reductions in gene flow among isolated populations. The European Tree Frog, Hyla 
arborea, has suffered from dramatic population declines in the last decades and is categorized as Vulnerable 
to Critically Endangered in its north-western distribution range. In Lower Saxony (Germany), the current 
distribution of the tree frog is fragmented. In this study, we aimed to assess the population structure, genetic 
diversity, gene flow, and migration rates in order to define the units for conservation management. Across 
a distribution area of 250 km’, frogs were sampled at 14 localities and genotyped at seven microsatellite 
loci, and the mtDNA cytochrome b gene was sequenced for a subsample. Whereas microsatellite pairwise 
D, „and F,, values showed genetic differentiation among nearly all sampled populations, Bayesian analyses 
assigned the 14 localities to two distinct genetic clusters including seven subclusters. Together with a 
slight correlation between geographic and genetic distance, the population structure indicates ongoing 
fragmentation. The cytochrome b haplotype distribution does not indicate divergence into mtlineages, but 
highlights the former connection of populations along the river Elbe. The results of this study suggest that the 
intense anthropogenic pressures in this area over the last decades have had negative genetic consequences 
for this species. The fragmented population structure calls for reconnection of the isolated occurrences by 
the implementation of conservation measures. 
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When reconnection of the habitats of endangered 


Introduction 


Genetic diversity and connectivity mediated by migrating 
individuals between populations are critical for the 
maintenance of many threatened species and can be 
evaluated by population and landscape genetic analyses 
(Shaffer et al. 2015). Loss of connectivity disrupts gene 
flow between formerly connected habitats and leads to 
the isolation of populations. Isolation in turn imposes a 
more rapid erosion of genetic diversity, exacerbating the 
effects of genetic drift and inbreeding on local gene pools 
(Andersen et al. 2004; Crnokrak and Roff 1999; Hedrick 
and Kalinowski 2000; Luquet et al. 2011). 


species is necessary, it 1s essential to determine the 
genetic structures and migration patterns for effective 
conservation management. This information can be used 
to delineate conservation units (e.g., Palsboll et al. 2007), 
even though the concepts that are applied to define them 
are somewhat uneven among studies and taxa (Shaffer 
et al. 2015). While Evolutionary Significant Units 
(ESUs) are used to delineate entities which possess a 
long (evolutionary) history (Crandall et al. 2000; Moritz 
1994), management efforts are often restricted to a more 
recent and regional space. In such cases, population 
boundaries need to be identified among which gene 
flow 1s limited. To achieve this, both mitochondrial 
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and nuclear markers are informative. Mitochondrial 
DNA has been widely used to analyze the phylogenetic 
relationships of amphibian populations (Dufresnes et 
al. 2013; Stock et al. 2012), while nuclear markers like 
microsatellites are well suited for detecting fine-scale 
structuring of populations and recent loss of genetic 
variation (Selkoe and Toonen 2006). Population genetic 
approaches such as Bayesian assignment tests use 
population allele frequencies to group individuals into 
genetic clusters. Together with information on genetic 
divergence between genetic clusters, this approach can 
be used to denote conservation units (Olsen et al. 2014; 
Rowe and Beebee 2007). 

In Europe, habitat loss, fragmentation, and degradation 
— mostly due to anthropogenic pressure (Cushman 2006; 
Pimm and Raven 2000) — are the most significant threats 
to endangered wildlife populations (Fahrig and Merriam 
1994; Sih et al. 2000; Stuart et al. 2004). Amphibian 
populations are especially vulnerable to fragmentation 
and loss of genetic variation due to their low dispersal 
capabilities (as reviewed in Smith and Green 2006). For 
safeguarding vulnerable species of this most endangered 
vertebrate group (Stuart et al. 2004), it is necessary to 
counteract genetic depletion by maintaining the exchange 
of individuals among populations. 

The European Tree Frog has shown long-term decline 
in much of its Western European distribution, mainly 
caused by habitat fragmentation (Andersen et al. 2004; 
Dubey et al. 2009; Krug and Pröhl 2013). The highest 
genetic diversity of this species has accumulated in 
South-eastern Europe, where it survived in refugia during 
glaciations. After late-Pleistocene diversification on the 
Balkan Peninsula, one of several major genetic groups 
recolonized North and Western Europe. Postglacial 
expansions resulted in decreasing genetic diversity across 
the range and therefore increased the vulnerability of 
populations towards North-Western Europe (Dufresnes 
et al. 2013; Stock et al. 2012). Indeed, the tree frog is 
not categorized as threatened in South-Eastern Europe, 
while it is reported to have declined and is now classified 
as Vulnerable to Critically Endangered in different 
areas in the north-west (see review in Dufresnes et al. 
2013, Table S1). In Lower Saxony in Germany, where 
the current distribution is patchy with some main 
occurrences in the lowlands (Fig. 1), the conservation 
status of the tree frog is Endangered (see the Red List at 
http://www.amphibienschutz.de). Although the species 
was widespread in the past, severe declines have been 
observed mainly in the second half of the last century 
(Manzke and Podloucky 1995). In some places, measures 
for conservation management have been successful 
(Brandt and Lüers 2017; Buschmann et al. 2006; Richter 
and Mugge 2012). 

For supporting further conservation activities, 
analyses of the genetic structure are required for 
assessing the genetic clusters as a way to delineate the 
units for conservation management (Rowe and Beebee 
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2007). So far, several studies have measured the genetic 
structure and diversity in more or less fragmented 
metapopulation systems (Andersen et al. 2004; Angelone 
and Holderegger 2009; Arens et al. 2006; Dubey et al. 
2009; Edenhamn et al. 2000; Krug and Pröhl 2013). The 
aim of this study was to perform a conservation genetic 
survey of the European Tree Frog across its distribution in 
Lower Saxony and adjacent areas. The specific intention 
was to assess significant genetic differentiation in order 
to define those conservation units among which dispersal 
is restricted. The obtained information was then used to 
identify population management goals and to provide 
specific recommendations about conservation priorities 
to ensure the long-term survival of the tree frog in this 
region. 

In this study, we tested hypotheses regarding 
population genetic structure, differentiation, and diversity 
in the Endangered European Tree Frog by analyzing 
mitochondrial sequence and nuclear microsatellite data 
with a series of statistical techniques. First, we tested 
for the existence of diverged genetic lineages. We 
further predicted that past population expansion and 
recent habitat fragmentation 1) reduced the migration 
among localities; 2) reduced genetic diversity as well as 
genetic population size within localities; and 3) resulted 
in significant genetic structure among the remaining 
tree frog localities. Therefore, 4) we expected a small 
to moderate effect of geographic distance on genetic 
differentiation as a result of the ongoing population 
disconnection. This work reveals ongoing population 
fragmentation with moderate genetic diversity for the 
populations of the European Tree Frog in Lower Saxony. 


Materials and Methods 
Sample Collection and Preparation 


Fourteen sites were sampled across the tree frog 
distribution in Lower Saxony and adjacent distributions 
in North Rhine Westphalia and Saxony-Anhalt, all in 
Germany. We chose one sample site within each main 
occurrence of the tree frog in this region (Fig. 1). In the 
occurrence near Hannover, however, we sampled four 
sites: two in the west of Hannover (KZ, KO, see Table 
1 for site definitions) and two in the east of Hannover 
(KH, BH) for a comparison of smaller scaled spatial 
distances. In total, 237 individuals were sampled with 
5-22 individuals per sample site (Table 1). Genetic 
material was collected from the tips of tadpole tails and 
buccal swabs of adult frogs. The adults were collected 
from the choruses during the breeding seasons in spring 
2007 and 2008. Tadpoles were sampled in summer 2007 
at three localities. In this year, the climatic conditions for 
breeding where unfavorable and adult catch rates were 
low at these sites. To avoid bias in the results from tadpole 
samples representing offspring from only one breeding 
pair, tadpoles were sampled in different breeding ponds. 
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Fig. 1. Current distribution of the European Tree Frog, and the distribution of cyt b haplotypes in Lower Saxony and adjacent areas 
on the basis of TK25-quadrants (grey squares) during 1994-2010 in Lower Saxony (NLWKN 2011), 1993-2006 in North Rhine 
Westphalia (LANUV 2011), and 1990-2000 in Saxony Anhalt (Meyer et al. 2004). Dashed lines denote state borders, dots denote 
sample sites. Inset in upper left corner: Haplotype network of 11 distinct haplotypes of cyt b of H. arborea (901 bp) in Lower 
Saxony and adjacent areas. Each haplotype is represented by one circle and color. The size of a circle corresponds to the haplotype 
frequency. Lines between haplotypes denote mutational steps between sequences. 


The genetic diversity for these localities was similar to 2.2. The genotyping results can be found in the file that 
other locations, suggesting that relatedness among samples ^ accompanies this article (Supplementary file 1). 
did not bias the results (Table 1). DNA from the tail clips Because earlier analyses (Stock et al. 2012) included 
was fixed in 99% ethanol and extracted using a proteinase only four samples of mt DNA from Germany, we also 
K digestion followed by a Phenol-Chloroform protocol ^ sequenced cytochrome b (cyt b) fragments of 901 bp for 
(Sambrook et al. 1989), and then stored at -20 °C. DNA 5-20 individuals from each sample site, excluding KO and 
was extracted from the buccal swabs with an Invisorb Spin BH from the Hannover population. The cyt b fragment 
Swab Kit (Invitek) following the manufacturer's protocol, was amplified via PCR using the primers MVZ 15-L (5'- 
and stored at -20 °C. Another study confirmed that buccal ^ GAACTAAT GGCCCA CACWWTACGNAA -3') and 
swabbing is a very efficient method for obtaining DNA of cyt b AR-H (TAWAAGGGTCTTCTACTGGTTG) from 
adequate quality for microsatellite amplification (Broquet Moritz et al. (1992) and Goebel et al. (1999). The PCR 
et al. 2007). reaction (25 ul) consisted of 20-100 ng DNA, 1 ul of 
A total of seven polymorphic microsatellite loci each primer (10 uM), 0.8 ul dNTPs (10 mM, SPRIME), 
(WHA1-9, WHA1-20, WHA1-25, WHA1-67, WHAI- 2.5 ul 10x advanced Buffer (SPRIME), 1.25 U Taq DNA 
103, WHA1-104, and WHA1-140) previously isolated Polymerase (SPRIME), and 17.45 ul H,O. The PCR 
by Arens et al. (2000) were amplified following the conditions were as follows: an initial denaturation at 94 °C 
author's protocol, except that the annealing temperature for 3 min; 35 cycles at94°C for 45 s, annealing temperature 
for WHA1-20 was changed to 64.6 °C. The PCR of50 °C for 45 s, and extension at 65 °C for 1 min. The 
products were genotyped using the capillary sequencer PCR products were sent to the Macrogen Company 
MegaBace 1000 (Amersham Bioscience). Allele scoring (Seoul, South Korea) for purification and sequencing with 
was performed using the software Genetic Profiler v. | an ABI3730XL genetic analyzer (Applied Biosystems). 
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Table 1. Overview of data from the various sample sites. *: Samples from adult frogs, ': samples from tadpoles, H; observed 
heterozygosity, H,: expected heterozygosity, SD: standard deviation, F: inbreeding coefficient, with bold values for significant 
differences after 1,000 permutations, R: mean allelic richness over all loci, A: haplotype diversity, x: nucleotide diversity, N: number 


of sampled individuals, N,: mean number of alleles over all loci. 


ID Sample site mean H, + SD mean H, + SD p IN R h x [96] N 
QU Quakenbrück ° 0.786 + 0.248 0.741 + 0.143 -0.070 4.00 — 0.00 0.00 5 
WK Westerkappeln? 0.661 + 0.187 0.579 + 0.158 -0.154 3.43 — 0.54 0.06 8 
EK Espelkamp: 0.796 + 0.169 0.754 + 0.087 -0.058 5.20 5.18 0.00 0.00 12. 
KZ Kananohe Zentrum ° 0.667 + 0.169 0.666 + 0.134 -0.001 5.14 4.65 0.00 0.00 20 
KO Kananohe Ost? 0.701 + 0.208 0.684 + 0.089 -0.027 4.57 4.53 _ — 11 
KH Kolshorn è 0.754 + 0.131 0.713 + 0.094 -0.059 6.29 337 0.53 0.13 20 
BH Beinhorn ° 0.693 + 0.089 0.693 + 0.093 -0.001 3.57 4.84 — — 20 
BA Bassum t 0.771 + 0.099 0.748 + 0.092 -0.032 5.43 4.98 0.41 0.05 20 
RU Ruschwedel * 0.731 + 0.064 0.721 + 0.050 -0.015 5.00 4.50 0.68 0.14 18 
WG Wolfsburg-Gifhorn? 0.790 + 0.158 0.799 + 0.080 -0.011 T1 6.65 0.61 0.08 20 
ST Strothe ^! C7850 153 0.708 + 0.126 -0.039 6.29 5.34 0.19 0.02 21 
AN Amt Neuhaus ° 0.708 + 0.111 0.750 + 0.090 0.057 6.43 5.50 0.66 0.15 22 
SW Salzwedel t 0.600 + 0.227 0.687 + 0.181 0.130 6.00 5.01 0.57 0.11 20 
PW Pevestorfer Wiesen? 0.793 + 0.110 0.764 + 0.091 -0.039 6.43 338 0.42 0.05 20 


Statistical Analysis 
Analysis of mtDNA 


Both directions of the cyt b sequences were assembled 
using the computer software SeqMan™ II (DNASTAR, 
Inc., Konstanz, Germany). Multiple sequence 
alignments were performed in MEGA 4 (Tamura et al. 
2007) using the Muscle algorithm (Edgar 2004), and 
all variable sites were confirmed by visual inspection 
of the chromatograms. The EMBL-EBI sequence 
analytical tool (Madeira et al. 2022) was used to 
convert the sequences to the corresponding amino acid 
sequences in order to assure that nuclear copies were 
not sequenced. The program MEGA was applied to 
calculate p-distances between sample sites (Tamura 
et al. 2004). Haplotype diversity (A) and nucleotide 
diversity (m) (Nei 1987) were determined with 
ARLEQUIN ver. 3.11 (Excoffier et al. 2005). A haplotype 
network of the cyt b data set was constructed via the 
Statistical parsimony analysis of the program TCS 
1.21 (Clement et al. 2000) using the default settings. 


Analysis of Microsatellites 


Microsatellite data were checked for null alleles, 
stuttering, and allelic dropout using Micro-CHECKER 
(Van Oosterhout et al. 2004). The program Fstar v. 
2.9.3 (Goudet 1995) was used to test for genotypic 
disequilibrium of all pairs of loci in each sample and to 
calculate average allelic richness per population. For the 
calculation of average allelic richness, sample sites with 
less than ten individuals (QU and WK) were excluded. 
For each sample site and locus, the observed and 
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expected heterozygosity (Nei 1987) and deviation from 
Hardy-Weinberg equilibrium (HWE) (Guo and Thompson 
1992) were determined with ARLEQUIN ver. 3.11 (Excoffier 
et al. 2005). Genepop ver. 4.1 (Rousset 2008) was used to 
test for a global deviation from HWE in each sample site. 
The inbreeding coefficient F, per sample site (Weir and 
Cockerham 1984) was calculated using Genetix ver. 4.05 
(Belkhir et al. 2004) and the significance was tested with a 
permutations test (1,000 permutations). 

Genetic differentiation between the sample sites was 
calculated as global F „and pairwise F „values (Weir and 
Cockerham 1984) in ARLEQuIN (Excoffier et al. 2005). In 
addition, pairwise D, (Jost 2008), a substitute measure 
of genetic differentiation, was calculated using the R 
package DEMEtics (Gerlach et al. 2010). Significance 
was calculated by 10,000 bootstraps. 

The data were also tested for Isolation By Distance 
in sampled populations (IBD; Storfer et al. 2010; Wright 
1943). IBD occurs when gene flow occurs but declines 
with increasing distances between pairs of populations, 
and is typical for the genetic population structure of 
many animal species (Hitchings and Beebee 1997; Spear 
et al. 2005; Vergara et al. 2015). To test for IBD, a Mantel 
test for correlation between pairwise genetic distances 
(F,, and D „) and pairwise geographic distances was 
conducted, implemented in IBDWS 3.23 (Jensen et al. 
2005). As proposed by Rousset (1997) for populations 
in two-dimensional habitats, geographical distance was 
log-transformed and genetic distance was expressed 
asel Jd rand. «(1D y Sıgnılcance for 
r 2 0 was assessed via 10,000 bootstraps. The linear 
geographic distances among sample sites were calculated 
in ArcView GIS 3.3 using the Distance Matrix extension 
(Jenness 2005). 
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Even though the Mantel test is widely used in 
landscape genetic studies, an evaluation of different 
methods revealed that Mantel tests exhibit high type- 
1 error rates (Balkenhol et al. 2009). Those authors 
recommended applying a combination of statistical 
methods to avoid inaccurate conclusions derived from 
only one method. Therefore, two additional hierarchical 
Bayesian methods, GESTE (Foll and Gaggiotti 2006) and 
BIMr (Faubet and Gaggiotti 2008), were applied here for 
evaluating the effect of distance by means of generalized 
linear models. Both BIMr and GESTE perform well for 
moderate samples sizes and limited numbers of loci, as in 
our study (Balkenhol et al. 2009). 

GESTE estimates the genetic distance (7... values) 
for each local population pair from multilocus genotypes 
and correlates them to environmental factors. Posterior 
probabilities associated with each factor allow the 
identification of factors with the highest effect on genetic 
structure. The regression coefficient estimate (Alpha) 
indicates whether a factor reduces or enhances genetic 
differentiation. The estimation of model parameters is 
performed by using a combination of Markov Chain 
Monte Carlo (MCMC) and Reversible-Jump MCMC 
(RJMCMC) (Green 1995). As environmental factors, 
we included latitude (G1) and longitude (G2; geographic 
coordinates in GK3 format) as approximations of the 
effect of distance among population pairs. 

The software BIMr 1.1 estimates contemporary gene 
flow and assesses the influence of genetic distance on gene 
flow. This program quantifies the gametic disequilibrium 
from multilocus genotypes (here, microsatellite alleles) 
generated by the progeny of recent migrants to calculate 
the proportion of the population that immigrated during 
the last generation (Faubet and Gaggiotti 2008). Five 
replicates (= runs) were run with a total of 1,020,000 
iterations (burn-in: 1,000,000, sample size: 20,000) and 
a thinning interval of 50 iterations. For each replicate, 
first 20 short pilot MCMC runs of 1,000 iterations 
were conducted, and the run with the lowest Bayesian 
deviance (D and the highest posterior probability 
was selected to extract the parameter estimates (Faubet 
et al. 2007; Faubet and Gaggiotti 2008). Two models 
were calculated: model 0 did not include environmental 
factors and model 1 included factor G1 which is the 
geographic distance between pairs of populations. As 
an alternative, the BAvEsAss software (Rannala 2007; 
http://www.rannala.org/software/) was also used to 
infer contemporary migrations rates. The software 
was run with 10,000,000 iterations (1), a burn-in (b) of 
1,000,000 repetitions, and the interval between samples 
(n) was set to 1,000. The default values were used for all 
other parameters at first. Then, we adjusted the mixing 
parameters for migration rate (m), allele frequencies 
(a) and inbreeding coefficients (f) to maintain. their 
acceptance rates between 20% and 60% as recommended 
in the Manual. 
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To infer genetic clusters, individual assignments to 
populations were conducted by means of a combination 
of non-spatial and spatial Bayesian algorithms with 
STRUCTURE version 2.3.3 (Pritchard et al. 2000) and TESS 
version 2.3 (Chen et al. 2007; Frangois et al. 2006). 
Simulation data suggested the combination of TESS and 
STRUCTURE as a reliable approach for deducing the spatial 
population structure (Chen et al. 2007), outperforming 
other Bayesian clustering programs. All STRUCTURE runs 
used 500,000 iterations after a burn-in period of 100,000. 
An admixture ancestry model and correlated allele 
frequencies were used between populations. STRUCTURE 
was run both without and with information about the 
sampling location (prior population information) and the 
results were compared as recommended by Pritchard et al. 
(2000; see also Dufresnes et al. 2013; Olsen et al. 2014). 
Hierarchical analyses were performed by repeating the 
STRUCTURE runs with each of the major clusters. Twenty 
runs were conducted for each K. The range of possible 
Ks tested spanned from 1 to 14, according to the number 
of sampled breeding sites. The average log likelihood 
Pr(X|K) (given by the estimated Ln Prob of data = Ln 
P(D) in the software result output, see Table 4) was 
calculated for each K across all runs. Since detecting 
the true number of K is not always straightforward, we 
included the AK statistics proposed by Evanno et al. 
(2005), using STRUCTURE HARVESTER v.0.6.8 (Earl and von 
Holdt 2012). 

TESS uses a Bayesian method to detect population 
structure, but it considers the spatial information 
(geographical coordinates) of the individuals. After 
assessing the preliminary runs as recommended in the 
software manual, the maximum number of allowed 
genetic clusters (Ka) was varied from 2 to 10. One 
hundred independent runs for each K „a were conducted 
under the admixture model, with 50,000 sweeps and a 
burn-in period of 10,000 sweeps for each run. 

Finally, the NEESTIMATOR v2 software was used 
to estimate the contemporary effective population 
sizes at all sample sites (Do et al. 2014). Three single 
sample estimators were implemented: the linkage 
disequilibrium method, the heterozygote-excess method, 
and the molecular coancestry method. The lowest allele 
frequencies (P...) were set to 0.05, 0.02, 0.01, and 0+. 


Crit 


Results 
Genetic Diversity 


Genetic diversity was estimated for each population based 
on the microsatellite alleles and the cyt b haplotypes (Table 
1). While genetic diversity indices based on microsatellites 
(H,, H,, N,, and R) are moderate to high across the range, 
they are always highest in WG (except for H ); while the 
indices based on cyt b (h, x) tend to increase from west 
(WK/QU) to east (PW) (Table 1; Fig. 1). 
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Table 2. Pairwise D, values (lower matrix) and pairwise F „ values (upper matrix) between sample sites; ns = not significant. See 


Table 1 for sample site acronym definitions. 


QU WK EK KZ KO KH 

QU 0 0.120 0.066 0.111 0.113 0.110 
WK 0.208 0 0.114 0.064 0.107 0.117 
EK 0.195 0.286 0 0.041 0.031% 0.069 
KZ 0.284 0.145 0.102 0 0.0025 0.090 
KO 0.324 0.243 0.097 0.001” 0 0.092 
KH 0.412 0.301 0.222 0.279 0.290 0 
BH 0.387 0.372 0.230 0.284 0.311 0.059 
BA 0.189 0.230 0.220 0.199 0.189 0:172 
RU 0.424 0.449 0.285 0.382 0.386 0.173 
WG 0.139 0.309 0.296 0.230 0.306 0.237 
ST 0.299 0.208 0.329 0.281 0.308 0.175 
AN 0.265 0:315 0.230 0.248 0.263 0.245 
SW 0.238 0,394 0.289 0.233 0.228 0.368 
PW 0.213 0.365 0.278 0.278 0.322 0.350 


BH BA RU WG ST AN SW PW 
0.103 0.039 0.117 0.043 0.086 0.086 0.084 0.052 
0.154 0.097 0.191 0.107 0.083 0.119 0.162 0.147 
0.076 0.060 0.090 0.061 0.095 0.065 0.097 0.071 
0.099 0.066 0.139 0.081 0.092 0.081 0.095 0.105 
0.104 0.059 0.130 0.081 0.092 0.074 0.091 0.106 
0.025 0.053 0.062 0.064 0.056 0.074 0.110 0.094 

0 0.074 0.092 0.082 0.076 0.102 0.118 0.100 
0225 0 0.074 0.037 0.072 0.070 0.083 0.066 
0232 0.235 0 0.068 0.079 0.072 0.110 0.073 
0.283 0.142 0:252 0 0.056 0.071 0.086 0.061 
0.238 0.239 0.222 0.201 0 0.057 0.110 0.070 
0.316 0.242 0.238 0.290 0.211 0 0.097 0.065 
0.367 0:272 0.347 0.287 0.369 0.343 0 0.078 
0.337 0.253 0.267 0.249 0.226 0.241 0.253 0 


Mitochondrial Sequence Analysis 


The analysıs revealed 11 haplotypes of the cytochrome 
b fragment which differed by ten variable sites and 
nine parsimony informative sites (Fig. 1). There was no 
evidence for any diverged haplotype groups that would 
correspond to different genetic lineages. Most haplotypes 
were closely related but unique to one sample site, except 
for haplotypes Hy-1, Hy-2, and Hy-5. While Hy-1 (blue) 
and Hy-5 (red) showed a broad distribution over almost 
the complete sampling area, Haplotype Hy-2 (green) 
was restricted to five sample sites in the northeast (Fig. 
1). Eight haplotypes were found at only a single locality: 
Hy-3 (white) in AN, Hy-4 (orange) in BA, Hy-6 (light 
blue) in EK, Hy-7 (yellow) and HY-8 (brown) in KH, 
Hy-9 (dark blue) in RU, and Hy-10 (grey) and Hy-11 (dark 
grey) in WG. In WG and RU, four different haplotypes 
were detected, while in QU, EK, and KZ (all in the west of 
Hannover) only one haplotype was found. The p-distances 
among localities were low, varying between 0 and 0.4 
% (Supplementary Table SI). The GenBank accession 
numbers can be found in Supplementary Table S3. 


Microsatellite Analysis 


The seven microsatellite markers examined were 
polymorphic with seven to 16 alleles per locus. The 
analysis with MICRO-CHECKER uncovered signs of null 
alleles for locus WHA1-67 in sample site KO and for 
locus WHA1-140 in sample site SW. As null alleles for 
these two loci were found at only a single sample site, we 
did not adjust for null alleles. Furthermore, this analysis 
revealed no evidence for large allele dropout or scoring 
errors due to stuttering. 

Deviation from Hardy-Weinberg-Equilibrium (HWE) 
was found for WHA1-104, with a significant excess of 
heterozygotes in sample site KH. The global test for 


Amphib. Reptile Conserv. 


HWE over all loci in each population resulted in no 
significant deviation from HWE. Significance values 
for the inbreeding coefficient F, were obtained for the 
sample sites SW (Fp = 0.130) and WK (F,, = -0.154, 
Table 1). No linkage (genetic) disequilibrium was found 
between any pair of loci. 

The global F „value across all localities was 0.083 and 
highly significant (P > 0.0001). Genetic differentiation 
calculated as pairwise D, and pairwise F „values were 
significant in all cases except between the two sample 
sites in the West of Hannover (KZ and KO), as well as EK 
and KO regarding the F... values (Table 2). The Mantel 
test for IBD showed a significant but low correlation 
between the genetic and geographic distances (Fig. 2; 
Di xr 02852 0.0197 eer 02952 — 00145 See 
also Supplementary Table S2), indicating that genetic 
differentiation is only partially explained by geographic 
distances among the sites. 

GESTE calculated five different models (Table 3). The 
probability of a model was not improved by including 
either latitude (G1) or longitude (G2) without interaction. 
The model with the highest posterior probability was model 
4, which included the constant, latitude, and longitude as 
well as their interaction. The Alpha values were low for the 
effects of both factors, while the Alpha value of the inter- 
action indicates a significant effect on genetic differentiation. 

All five replicates of the BIMr analysis showed a 
D of 0.0. The highest posterior probability for the null 
model was 0.79 (run 1), and the lowest was 0.55 (run 
4). The posterior probabilities for model 1 (including 
G1) were lower than the posterior probabilities of the 
null model (Table 4); i.e., the geographic distance did 
not seem to affect recent gene flow or migration among 
sample sites. Mean migration rates were extremely low 
and varied from 2.88e1? to 1.11e?; while the highest 
mean migration rate was observed among PW and 
WG, and the lowest was among AN and QU. Also, the 
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migration rates between two close population pairs in the 
surroundings of Hannover were very low. For example, 
the mean migration rates between KH to BH and vice 
versa were only 1.66e'? and 2.34e'°, and those between 
KZ and KO were 1.92€? and 2.24€? even though the 
distances between these pairs were only 3.94 km and 2.6 
km, respectively. The analysis with BAvEsAss provided 
very similar results for the runs with default and adjusted 
parameters. The migration rates were higher than those 
calculated with BIMr, most of which fluctuated around 
0.01. Interestingly, the migration rates between EK and 
KZ, KO and KZ, as well as KH and BH were considerably 
higher (0.14, 0.16, and 0.20, respectively) and more 
consistent with the genetic population structure (Fig. 3) 
as well as the above-mentioned geographic distances. 
Bayesian assignments conducted by STRUCTURE 
suggested two major groups (K = 2), separating the 
western/central populations (WK-BA-EK-KZ-KO-BH- 
KH) and the northern/eastern populations (RU-WG-ST- 
AN-PW-SW, Figs. 3A and 4, see also Supplementary 
Fig. SIA). The population QU appeared to be admixed, 
but was unambiguously assigned to the western group 


by TESS (see below). Within each cluster, a fine 
substructuring could be detected, and the two approaches 
(with/without prior population information) provided 
slightly different clustering solutions. In those runs where 
the sampling location was used as prior information, the 
western/central cluster was split along Hannover with 
some admixed populations in the western part of the 
region (K =3, Fig. 3B). The STRUCTURE runs without prior 
population information supported the existence of two 
genetic subclusters within this region (Supplementary 
Fig. SIB). In the northern/eastern cluster, both approaches 
provided similar results (Fig. 3C and Supplementary Fig. 
S1C). RU was differentiated from WG and ST, which in 
turn differed from the populations lying in the eastern 
part of the study area in Saxony Anhalt. More detailed 
results of log likelihood Pr(X|K) values and AK statistics 
are provided in the Supplementary figures (Figs. S2 and 
S3, respectively). 

The TESS analysis supported the two major genetic 
groups separated into a western/central cluster and a 
northern/eastern cluster (Fig. 5). Increasing the number 
of K „x resulted in only a slight decrease in the Deviance 


Table 3. Posterior probabilities for five possible models calculated with GESTE explaining the genetic differentiation of European Tree 
Frogs as a function of the environmental factors latitude (G1) and longitude (G2). Constant is the intercept of the regression model. The 
regression coefficients (Alpha) for different environmental factors used in the models are given in the right side of the table. 


Model Factors included Posterior probability Factors Regression coefficients 
ModelO Constant 0.19 Constant Alpha 0 -3.69 
Model 1 Constant, G1 0.19 G1 Alpha 1 0.06 
Model2 Constant, G2 0.18 G2 Alpha 2 -0.07 
Model 3 Constant, G1, G2 0.17 G1*G2 Alpha 3 -2.14 
Model 4 Constant, G1, G2, G1*G2 0.26 
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Fig. 3. Estimation of the number of Hyla arborea populations using the program STRUCTURE ver 2.3.1 (Pritchard et al. 2000) for the 
admixture model with prior population information; QU, WK, EK, etc. = sample sites, separated by fine black lines. Each individual 
is represented by a single vertical line broken into colored segments, with lengths proportional to the corresponding clusters. (A) 
Plot for X = 2 in the analysis of the entire data set, (B) plot for X = 3, and (C) K = 4 for hierarchical analysis on each of the two 


main clusters. 


Information Criterion (DIC) of the models while the 
resulting population structure was not consistent among 
runs for each K pax or in comparison with the STRUCTURE 
results (data not shown). STRUCTURE as well as TESS 
revealed that frogs from WG in the eastern group (Fig. 
3A) are of admixed origin from both genetic groups as 
the assignment results from the two were similar. 


Effective Population Sizes 


Mean effective population sizes (N.) varied among the 
sample sites and statistical methods applied. The results 
also differed between Po, = 0.05 and the other P... but 
were the same for P „= 0.02, 0.01, and 0+ (Table 5). The 
large confidence intervals indicate that the results might 
not be very reliable. However, most calculated N, values 
were small (100 individuals or less). Only three sample 
sites (KO, AN, and PW) showed consistently high values 
for N, (500 or higher, or infinite) for most methods. 


Discussion 


The analyses presented here provide valuable information 
for the conservation management of the Endangered 
European Tree Frog species, Hyla arborea, which suffers 
from population isolation in its northern distribution 
range. In Lower Saxony, a weak correlation between 
genetic and geographic distances suggests a low level of 


recent gene flow among localities, and further analyses 
indicate a lack of current migration at least during the last 
generation. Two major genetic clusters, one in the east 
and one in the west, were found with some admixture 
in a central population. Both main clusters were further 
subdivided into several distinct regional clusters. The 
substantial population structure, verified by significant 
genetic distances among localities, suggests that the 
populations are currently isolated to a large extent. 
Consequently, conservation management is needed to 
ensure the long-term persistence of this species in Lower 
Saxony with suitable effective population sizes and high 
levels of genetic diversity that are necessary to counteract 
the reductions in fitness and adaptive potential (Andersen 
et al. 2004; Frankham 2005; Allentoft and O'Brien 2010; 
Angelone 2010). 


Isolation by Geographic and Genetic Distances 


In addition to the Mantel test, the landscape genetic 
analysis in GESTE provided some insight into the role 
of geographic distance on genetic differentiation. While 
the models including latitude and longitude alone did 
not offer a better explanation for genetic differentiation 
than the null model, the most complex model including 
the interactions between latitude and longitude did. We 
interpret this as the effect of the geographic distance 


Table 4. Results of BIMr analysis (means of posterior probability and Alpha for run 1) for estimating migration rates among 
localities with European Tree Frogs in Lower Saxony. The factor G1 is the geographic distance. Alpha 0 and Alpha 1 represent 
estimates of the constant term and factor G1, respectively. 


Model Factor Posterior probability Alpha 0 Alpha 1 
Model 0 0.79 1.64 
Model 1 G1 = geographic distance 0.21 1.48 0.14 
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Fig. 4. Mean values of estimated Ln probability of data (LnPD) for each K (a) and delta K (b) when prior population information 


was implemented. 


on genetic dissimilarity between localities. Overall, 
however, the posterior probability and alpha values 
illustrate that the distance effect is not very large. 
Habitat fragmentation might play a more important role 
in shaping the genetic structure of the tree frogs in this 
area. According to Podloucky and Fischer (2013), habitat 
fragmentation in Lower Saxony is mainly caused by the 
loss of summer habitat, breeding ponds, and corridors 
suitable for migration. This is in accordance with the 
disconnected distribution (Fig. 1) and significant genetic 
distances (F...and D, ) among most tree frog populations. 

A significant population structure as a result of 
limited dispersal between isolated populations is 
typically accompanied by a slight to moderate effect of 
isolation by distance. Isolation by distance was detected 
in some European amphibian species (Rana dalmatina, 


Fig. 5. Map of geographic-genetic cluster membership for K pax 
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Sarasola-Puente et al. 2012; Bombina variegata: 
Weihmann et al. 2009; Hantzschmann et al. 2020), but 
not others (Bufo calamita: Allentoft et al. 2009; Bombina 
bombina: Dolgener et al. 2012). For the tree frogs in 
Lower Saxony, the small positive correlation (r ~ 0.28) 
between genetic and geographic distances suggests 
a very low level of recent gene flow among localities. 
This finding is in accordance with earlier studies on 
tree frogs, which reported small to moderate correlation 
coefficients between both distances and significant 
population structure (Andersen et al. 2004; Angelone and 
Holderegger 2009; Arens et al. 2006). In all these studies, 
limited gene flow was explained by habitat fragmentation, 
particularly the loss of breeding ponds. In contrast, 
in those frog species which occur in more continuous 
habitats or that have higher dispersal capacities, the 


— 2 as inferred by TESS. 
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Table 5. Results for effective populations sizes (NV) calculated with the N, Estimator for three different methods and different P... 
values. N, values were the same for P... = 0.02, 0.01, and 0+. Only O+ is used in the Molecular Coancestry method. Mean values 


(N.) and Confidence Intervals (CIs) are given for each sample site, method, and P... option. 
The three methods 
Linkage disequilibrium Heterozygote excess Molecular 
coancestry 
Sample site p 0.05 0.02/0.01/0+ 0.05 0.02/0.01/0+ 0+ 
QU N, 35:7 357 7.8 7.8 9.5 
CIs 1.4 - inf. 1.4 - inf. 3.7 — inf. 3.7 — inf. SW —y52 
WK N, 1.3 1.3 Sal 5.1 2.7 
CIs 0.7-2.7 0.7-2.7 2.5 — inf. 2.5 - inf. 2.0-3.5 
EK N, 16.3 33:5 9.9 9.9 8.2 
CIs 5.6 — inf. 8.7 — inf. 4.5 — inf. 4.5 — inf. 2.7 — 16.9 
KZ N, 46.3 28.2 34.7 34.7 5.6 
CIs 15.8 — inf. 12.1 2 3192 7.5 — inf. 7.5 — inf. 2.3— 10.5 
KO N, inf. 611.3 231.5 2315 inf. 
CIs 11.4 — inf. 10.4 - inf. 4.3 — inf. 4.3 — inf. inf. — inf. 
KH N, 319.3 159.7 inf. 44.3 37 
CIs 30.5 — inf. 30.2 - inf. 5.8 — inf. 6.0 — inf. 3.2— 8.9 
BH N, 83.5 TIS 73157. inf. inf. 
CIs 21.0 — inf. 21.3 - inf. 7.2 — inf. 7.7 — inf. inf. — inf. 
BA N, 1772 188.0 322.0 322.0 10.8 
CIs 26.5 — inf. 26.8. inf. 6.5 — inf. 6.5 — inf. 1.8 —27.7 
RU N, 19.3 21.5 67.7, 81.9 inf. 
CIs 8.4 — 106.4 10.8 — inf. 6.8 — inf. 7.1 — inf. inf. 
WG N, 311 65.5 inf. inf. 9.1 
CIs 16.0— 113.5 26.0 — inf. 8.1 — inf. 9.3 — inf. 3.8— 16.6 
ST N, 34.7 TES 39.7 51.0 inf. 
CIs 157-3157 25.4 - inf. 6.4 - inf. 72- inf. inf. — inf. 
AN N, inf. inf. inf. inf. inf. 
CIs 41.7 — inf. 74.0 — inf. 30.9 — inf. 57.8 — inf. inf. — inf. 
SW N, 37.3 36.0 inf. inf. 43.4 
CIs 15.2 - inf. 16.1 - inf. inf. — inf. 280 - inf. 0—217.7 
PW N, 631.9 inf. inf. inf. inf. 
CIs 36.2 — inf. 73.1 — inf. 6.1 — inf. 6.6 — inf. inf. — inf. 


correlation between geographic and genetic distance 1s 
higher (e.g., Próhl et al. 2006) or disappears in case of 
panmixia (Leblois et al. 2000). 

The BIMr analysis revealed that contemporary 
migration (1.e., during the last generation) was practically 
absent among the sample sites of tree frogs in Lower 
Saxony. This is incontrastto the migration rates calculated 
with BayesAss and the results of the STRUCTURE analysis 
(see below), which imply that individuals of mixed 
ancestry exist in the different genetic clusters. However, 
gene flow over the last few tree frog generations may 
have suffered due to the expanding habitat fragmentation, 
while some decades ago far more tree frog localities were 
reported for lower Saxony, thus, the connectivity among 
them was much better (Manzke and Podloucky 1995). 
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Therefore, the signals for gene flow are still apparent in 
the results of some analyses (STRUCTURE, BAYESAss) but 
not in those where the calculations are restricted to the 
most recent years (BIMr). 


Genetic Diversity 


The expected microsatellite heterozygosity (H ) has been 
measured in a number of previous population genetic 
studies of the European Tree Frog. Interestingly, H, was 
higher in the current study area (H : 0.60 — 0.79) than in 
most other northern areas where the genetic situation of 
the tree frog was investigated (e.g., H, values in Denmark: 
0.35—0.54, Andersen et al. 2004; Switzerland: 0.27—0.71, 
Angelone and Holderegger 2009; and the Netherlands: 
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0.39-0.59, Arens et al. 2006). As expected, the more 
peripheral populations, such as those in Denmark and 
the Netherlands, show lower genetic diversity values. In 
comparison, the mthaplotypes diversity (mean h=0.38, n 
= 14 populations) is lower in the Lower Saxony area than 
in the southern part of the distribution range (Greece, 
Albania, Croatia, Serbia, and Romania) were the average 
h amounts to 0.7 (n = 20 populations, calculated from 
Dufresnes et al. 2013, Table S2). From all cumulative 
data within the framework of this study, we can conclude 
that ın central and northern Europe, human induced 
fragmentation processes involving habitat destruction 
in a previously widely distributed frog species are 
contributing to the depletion of genetic diversity. 


Genetic Structure and Conservation Units 


The Bayesian cluster analyses conducted with STRUCTURE 
and TESS support the division of the tree frog populations 
into two major geographic-genetic clusters, one in the 
west and one in the east. One population in the south- 
east of the area (WG) shows admixture between the two 
groups and therefore a relatively high genetic diversity. 
Thus, itseems that both groups were previously connected 
by migrating animals when habitat fragmentation was 
less severe. This result is consistent with the moderate 
correlation between geographic distance and genetic 
distances. Both analyses also provide evidence for 
further fragmentation within both groups. In contrast, the 
haplotype network does not indicate any older, distinct 
(e.g., postglacial) lineages supporting the finding of 
Dufresnes et al. (2013), that only one evolutionary unit is 
present in this area. 

Microsatellite pairwise D, and F values all showed 
significant genetic differentiation except forthetwo closest 
sites (KZ and KO) in the West of Hannover. However, 
mthaplotype distribution and Bayesian analyses of the 
microsatellites suggest distinct relationships among the 
currently fragmented localities. In the Northeast of Lower 
Saxony, the distribution of the mthaplotypes indicates a 
former connection of the populations along the river Elbe. 
Interestingly, the easternmost occurrences at AN, PW, 
and SW in the current and former distributions display a 
relatively well-connected area, nonetheless the presently 
distinct genetic sub-clusters in the Bayesian analyses and 
pairwise F,, values are relatively high. The significant 
F, value found for SW indicates that the separation of 
this site may have resulted in inbreeding in an isolated 
population. Altogether, the available microsatellite data 
point to recent fragmentation of tree frog populations 
in this area, while similar mthaplotypes provide some 
evidence for a former connection. 

Interestingly, there ıs a significant genetic divergence 
between the sample sites KZ and KO in the West and 
the sample sites KH and BH in the East of Hannover. 
One possible explanation is that recently constructed 
motorways in combination with genetic drift contributed 
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to population differentiations, which are also apparent in 
the mthaplotype frequencies. Roads have been identified 
as barriers to gene flow in some other amphibians (Arens 
et al. 2007; Lesbarreres et al. 2006). One central question 
is whether these relatively young barriers (motorways 
expanding in the 1960s, and dense urban areas) are the 
only reason for the differentiation of these formerly 
linked tree frog localities (Manzke and Podloucky 1995). 
The low haplotype diversity in KZ, EK, and QU points to 
a loss of genetic diversity as a consequence of increased 
genetic drift in isolated occurrences. 

In summary, the genetic analyses point to a highly 
structured population, as was observed in other surveys 
of European amphibians (Dolgener et al. 2012; Rowe and 
Beebee 2007; Sarasola-Puente et al. 2012; Hantzmann et 
al. 2020). The risks of fragmentation include population 
reduction, loss of genetic diversity and declining fitness, 
and finally extinction (Hitchings and Beebee 1997; 
Cushman 2006). To reverse such negative processes, 
conservation management that takes the genetic 
population structure into account is important (Allentoft 
et al. 2009; Olsen et al. 2014). For this endangered and 
fragmented frog species, we suggest delineating the two 
major genetic clusters as conservation units; and then 
within those, intense reconnection efforts should be 
undertaken by creating suitable habitats for migrating 
frogs. Moreover, there is clear evidence of admixture in 
WG, and gene flow along this route could be reestablished 
between both clusters. 


Conservation Measures 


In our opinion, future conservation management should 
be directed towards two aims. The first aim is to maintain 
high genetic diversity in large and stable populations 
within each conservation unit. In this context, it has been 
argued that an effective breeding size (N) of at least 50 
animals is necessary to avoid inbreeding in the short 
term and that an N, of 500 is necessary to maintain the 
evolutionary potential that would allow adaptations to 
environmental changes and assure long term viability 
(Jamieson and Allendorf 2012). Most of our isolated 
sample sites (e.g., BA, Fig. 3B-C) or sub-clusters (e.g., 
KH-BH) do not reach these effective population sizes. 
For tree frogs, the ratio of effective breeding size (N) to 
census size (N) is ~ 0.5 (Broquet et al. 2009). Therefore, 
we recommend the monitoring of population sizes and 
maintaining population sizes of at least 100 breeding frogs 
in isolated populations, i.e., each sub-cluster, for short 
term conservation goals, but increasing the population 
sizes to 1,000 or more embedded in each of several meta- 
population systems in every conservation unit (see also 
Andersen etal. 2004; Frankham etal. 2014). In cases where 
populations within a conservation unit are genetically and 
geographically separated, genetic rescue can be attained 
by establishing corridors to stimulate dispersal (e.g., 
in AN-SW-PW). Dense networks of suitable spawning 
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ponds have been destroyed by habitat conversion but are 
of great importance for the maintenance of large tree frog 
populations and the connection of subpopulations. There 
are several reports that tree frogs not only respond well 
to new suitable water bodies, but also depend on them 
for migrations exceeding several km and often colonize 
them in subsequent breeding seasons (e.g., Angelone 
and Holderegger 2009; Brandt and Lüers 2017; Hansen 
2004; Schwartze 2002; Zumbach 2004). 

The second aim is to maintain overall genetic diversity 
among the genetic clusters within the species and to 
protect local co-adapted gene complexes (Savolainen et 
al. 2013). To achieve this second goal, we recommend 
the re-establishment of gene flow between genetic 
clusters where possible, but at a lower level than within 
them. This particularly applies to the western and eastern 
clusters between which (former) gene flow is evidenced 
by the STRUCTURE analysis (Fig. 3). Habitat reconnection 
between these areas would allow a few frogs to travel 
between the breeding ponds of different clusters, thereby 
refreshing genetic diversity and counterbalancing the loss 
of genetic diversity through drift, while diverse selection 
pressures would sustain local adaptation. Levels of genetic 
diversity inferred from neutral markers are not necessarily 
correlated with variation in locally adapted traits. In this 
context, more research is necessary to understand which 
traits are locally adapted and how their variation affects 
the fitness of a population. Very isolated and small 
populations might benefit from translocations, i.e., the 
introduction of individuals from other populations. In 
such cases, translocations should be restricted to within 
the conservation unit to avoid causing outbreeding 
depressions that have sometimes been observed between 
distantly related populations (Sagvik et al. 2005). 


Conclusions 


Populations of the European Tree Frog in Lower Saxony 
are highly fragmented geographically and genetically, and 
therefore endangered. We identified two major genetic 
clusters and recommend that they should be considered as 
local conservation units. Conservation efforts should entail 
areconnection of the populations within these conservation 
units, and to a lesser degree between them. Moreover, 
the maintenance of large and stable meta-populations 
within genetic sub-clusters (mostly consisting of isolated 
populations) needs to be achieved for long-term survival. 
For translocations of individuals to recovering very small 
and inbred populations or for reintroduction, we suggest 
a mixing of individuals from different populations within 
the same conservation unit to increase genetic diversity 
and enhance the adaptive capacity regarding changing 
environmental conditions. This study offers one example 
of how population genetic studies can help to delineate 
conservation units, and our recommendations might apply 
just as well to other endangered species where declines are 
connected to increasing habitat fragmentation. 
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Supplementary Fig. S1. Estimation of the number of Hyla arborea populations using the program STRUCTURE ver. 2.3.1 (Pritchard 
et al. 2000) without prior population information; QU, WK, EK, etc. = sample sites, separated by fine black lines. Each individual 
is represented by a single vertical line broken into K-colored segments, with lengths proportional to each of the K-inferred clusters. 
(A) Plot for K = 2 in the analysis of the entire data set, (B) plot for K = 4, and (C) K = 4 for hierarchical analysis on each of the two 
main clusters. 
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Supplementary Fig. S2. Mean values of estimated Ln probability of data (LnPD) for each K (A, C) and delta K (B, D) calculated 
from STRUCTURE runs with STRUCTURE HARVESTER (20 replicates per K) in those analyses where prior population information was 
implemented. (A-B), graphs for hierarchical analysis of the red cluster; (C-D), corresponding graphs for the green cluster (compare 
to Fig. 3 in the main text). 
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Supplementary Fig. S3. Mean values of estimated Ln probability of data (LnPD) for each X (A, C, E) and delta K (B, D, F) without 
prior population information. (A-B), results for the entire data set. (C-D), graphs for hierarchical analysis of the red cluster. (E-F), 
corresponding graphs for the green cluster. 
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Supplementary Table S1. Estimates of evolutionary divergence over cyt b sequence pairs between sample sites (p-distances). 


QU WK EK KZ KH BA RU WG ST AN SW PW 
QU 0 
WK 0.001 0 
EK 0.001 0.002 0 
KZ 0.000 0.001 0.001 0 
KH 0.003 0.002 0.004 0.003 0 
BA 0.001 0.001 0.003 0.001 0.002 0 
RU 0.001 0.001 0.002 0.001 0.003 0.002 0 
WG 0.001 0.001 0.002 0.001 0.002 0.001 0.001 0 
ST 0.001 0.001 0.002 0.001 0.002 0.000 0.001 0.001 0 
AN 0.001 0.001 0.002 0.001 0.003 0.001 0.002 0.001 0.001 0 
SW 0.001 0.001 0.002 0.001 0.003 0.001 0.002 0.001 0.001 0.001 0 
PW 0.001 0.001 0.002 0.001 0.002 0.001 0.001 0.001 0.000 0.001 0.001 0 
Supplementary Table S2. Geographic distances (km) among sample sites. 
QU WK EK KZ KO KH BH BA RU WG SI AN SW PW 
QU 0 24.89 59.36 118.51 121.07 138.45 138.18 59.60 145.74 186.98 194.94 225.27 214.32 244.29 
WK 24.89 0 55.80 121.60 124.18 140.04 140.47 742 165.05. 18968 20653, 23657 220551 723895 
EK 59.36 55.80 0 66.97 69.52 84.52 85.18 57.94 135.51 134.27 15858 187.71 168.42 203.12 
EZ 11551 12160 66.97 0 2.59 20.80 19.80 79:20: - T0770 68.55 98.15 125.14 101.84 138.06 
KO 121.07 124.18 69.52 259 0 18.53 17.33 81.32 107.77 65.97 9627 a 123.05 93397 135.79 
KH 138.45 140.04 84.52 20.80 18.53 0 3.94 99.8] 118.84 49.76 9508 117-07 88.33 126.83 
BH 138.18 140.47 85.18 19.80 17:33 3.94 0 98.08 115.08 49.22 89.62 114.01 86.17 124.28 
BA 59.60 77.42 57.94 79.20 81.32 99.81 98.08 0 87.65 142.67 137.64 167.94 161.85 188.28 
RU 145.74 165.05 135.51 107.70 107.77 118.84 115.08 87.65 0 134.01 84.08 108.84 12469 133.14 
WG 186.98 189.68 134.27 68.55 65.97 49.76 49.22 142.67 134.01 0 13.62 86.47 47.60 88.04 
ST 194.94 206.53 158.58 98.15 96.27 93.08 89.62 137.64 84.08 73.62 0 30.35 43.49 S17 
AN 22527723657 187.71 T2514- 23:05. 117.07- 11401. 167.94: 10834 86.47 30.35 0 42.06 24.53 
SW 21432 221.51 168.42 101.84 99.39 88.33 86.17 161.85 124.69 47.60 43.49 42.06 0 40.48 
PW 244.29 254.05 203.12 138.06 135.75 126.83 124.28 188.28 133.14 88.04 SUNT 24.53 40.48 0 


Supplementary Table S3. GenBank accession numbers for Hyla arborea CytB haplotypes Hy 1 to Hy 11. 


Sequence ID GenBank accession number 
BankIt2634361 Segl OP690610 
Bank1t2634361 Seq2 OP690611 
BankIt2634361 Seq3 OP690612 
Bank1t2634361 Seq4 OP690613 
Bank1t2634361 Seq5 OP690614 
Bank1t2634361 Seq6 OP690615 
BankIt2634361 Seq7 OP690616 
Bank1t2634361 Seq8 OP690617 
Bank1t2634361 Seq9 OP690618 
Bank1t2634361 Seq10 OP690619 
Bank1t2634361 Seq11 OP690620 
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